#########################################################
# Application
# Replication Rasmussen 2010
# Q1
#########################################################

library(readstata13)
source("cbq_binary_function.R")
    
# downloaded from http://journals.sagepub.com/doi/suppl/10.1177/1465116510388675
dat_r <- read.dta13("EUP388675-supplementary_material.dta")
# recode
dat_r$var15 <- as.numeric(as.character(dat_r$var15))
dat_r$var5 <- as.integer(dat_r$var5)
dat_r$var5[which(dat_r$var5==1)] <- 0
dat_r$var5[which(dat_r$var5==2)] <- 1
dat_r$var16 <- as.integer(dat_r$var16)
dat_r$var16[which(dat_r$var16==1)] <- 0
dat_r$var16[which(dat_r$var16==2)] <- 1
dat_r$var10 <- as.integer(dat_r$var10)
dat_r$var10 <- dat_r$var10 - 1 
dat_r$var15 <- dat_r$var15 - 1998
dat_r$var_int <- dat_r$var5 * dat_r$var6
dat_r$urgcy_yes <- ifelse(dat_r$var1=="Yes",1,0)
dat_r$urgcy_noleg <- ifelse(dat_r$var1=="No existing leg",1,0)
dat_r$var14 <- as.integer(dat_r$var14) - 1 
dat_r$var7 <- as.integer(dat_r$var7) - 1
dat_r$var17 <- as.integer(dat_r$var17) - 1
dat_r1 <- subset(dat_r,dat_r$var13==0)
dat_r1 <- subset(dat_r1,select = c(var16,
                                   var6,
                                   var3, 
                                   var5,
                                   var10,
                                   var15,
                                   urgcy_yes,
                                   urgcy_noleg,
                                   var4,
                                   var11,
                                   var14,
                                   var7,
                                   var17))
dat_r1 <- dat_r1[complete.cases(dat_r1),]
dat_r1$choice <- as.integer(dat_r1$var16)
# X,Y
Y <- dat_r1$choice
X <- as.matrix(subset(dat_r1,select = c(var6,var3, var5,var10,
    var15,urgcy_yes,urgcy_noleg,
    var4,var11,var14,
    var7,var17)))
N <- length(Y)
n_covariate <- dim(X)[2]
# fit the model
cbq_fit <- cbq(N = N,
               n_covariate = n_covariate,
               X = X,
               Y = Y,
               qtl = 1,
               nchain = 5,
               niter = 20000,
               thin = 20)
# save
save(cbq_fit,file="cbq_binary_q1.RData")